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PREDICTION  OF  THE  MOTION  OF  MISSILES 
ACTED  ON  BY  NON-LINEAR  FORCES  AND  MOMENTS 

ABSTRACT 

The  application  of  the  usual  techniques  of  non-linear  mechanics 
to  fourth  order  non-linear  systems  is  severely  handicapped  by 
algebraic  complexities.  It  is  shown  that  for  an  important  subset 
of  the  set  of  fourth  order  systems,  the  use  of  the  complex  variable 
allows  the  quick  derivation  of  the  required  results.  This  technique 
is  applied  in  some  detail  to  the  prediction  of  missile  yawing  motion. 
In  this  discussion  the  useful  concept  of  an  amplitude  plane  is  in¬ 
troduced.  Comparison  of  the  theory  with  the  results  of  exact  computa¬ 
tions  indicate  the  value  of  the  equivalent  liineanizdftton  approach. 

The  effects  of  gravity- induced  yaw  of  repose  and  of  small  aero¬ 
dynamic  asymmetries  on  the  general  non-linear  problem  Are  dismissed 
in  an  appendix. 


TABLE  OF  SYMBOLS 

A,B,C,D  complex  constants  in  Eq.  (2) 
a,b,c,d  real  constants  in  Eq.  (58) 


A  = 


V  - 


V  *2 


V 


wi 
D  = 
d 


real  coefficients  in  Eq.  (16) 
real  coefficients  in  Eq.  (16) 

real  coefficients  in  Eqs.  (l) 

,2 


(b  -  c )  +  4  ad 

diameter 


G  =  y< 


Srp 

H  = 

J.  = 


g. 


4-  ("d  '  ki2  V  -  Jg  +  & 

trajectory  component  of  the  gravitational  acceleration 
components  of  the  gravitational  acceleration  perpendicular  to 
the  missileifs  axis 


Pd" 

ra 

pd" 


KL  "  ^  +  k2  ^  “  iKMA^ 


Jg  = 


Ki  = 


m  Ki 


n 


z; 


,2k 


k  =  0 


2k 


i  =  H,  L,  M,  MA,  T 


PT 

Only  those  symbols  which  appear  in  the  body  of  this  report  are  listed 
here.  Symbols  which  are  introduced  in  the  appendices  appear  close  to 
their  definitions. 
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Ki  =  KiO  e 


-  aj_p 


i  =  1,  2 


amplitude  of  i  -  th  frequency 


% 

KL 

SlA 

K^p) 

*1 

*2 


moment  coefficient  due  to  cross  angular  velocity 
lift  force  coefficient  due  to  yaw 

moment  coefficient  due  to  yaw  (static  moment  coefficient) 
moment  coefficient  due  to  cross  acceleration 
Magnus  moment  coefficient 

parametric  function  for  damping  correction  to  first 
approximation 

axial  radius  of  gyration  in  calibers 
transverse  radius  of  gyration  in  calibers 
cosine  of  yaw  angle 


M  = 


m 


4  p 


m  = 

mi,m2 

P(x,y) 

P 

P 

Q(x,y) 


5 

m 


T  = 


ik2  “m 


mass 

complex  constants  in  Eq.  (5) 
function  in  Eq.  (58) 
independent  variable 

arclength  along  the  trajectory  in  calibers 
function  in  Eq.  (58) 


-2 


vki  y 

magnitude  of  velocity 

real  dependent  variables  in  Eqs .  (l) 
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0^  exponential  damping  coefficient  of  i-th  frequency 

(e2  +  ig,)d 

f  =  3  -  Jg  *■ 

u  6 


&  =  |/XX  (sine  of  yaw  angle) 


[>1  =  1_ 

L  el  2 it 


2it 

r 


Jo 


.2k 


K 

1  +  -5T-  COS  0 
K1 


a0 


I?*] 


e2 


6 


1 

2jt 


2« 

/“ 


.2k 


K, 


1  +  ~  cos  0  I  d0 
*2 


-] 


-'0 
arg  m^ 

complex  dependent  variable  in  Eq.  ( 3  );>or  complex  yaw 


1  “ld 

— ~  (gyroscopic  spin) 


5 

p 

a 


complex  dependent  variable  in  Eq.  (2) 
air  density 
arg  C 


0±  =  0^o  +  0^P  phase  angle  of  i-th  frequency  _ 


0  =  0X  - 

*2 

i(p) 

parametric  function  fdr  frequency  correction  of  first 

approximation 

O) 

arg  D 

®1 

axial  spin 

tilde  superscript  denotes  quantities 
epicyclic  first  approximation  of  the 

appearing  in  the 

non-linear  equation 
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1.  INTRODUCTION 


One  of  the  most  challenging  current  problems  in  exterior  ballistics 
is  the  prediction  of  missile  motion  for  missiles  flying  at  large  yaw 
angles.  Highberg  and  Zaroodny  have  treated  the  special  case  of 
circular  yawing  motion  while  the  case  of  almost  circular  motion  has  been 
rather  elegantly  analysed  by  Davis,  Follin,  and  BlitzerA  The  problem 
of  more  general  motion  for  cubic  static  and  Magnus  moments  has  been 
attacked  by  LietmamA  but  this  treatment  is  obscured  by  lengthy  algebraic 

manipulations  and  does  not  give  its  results  in  a  convenient  form  for  the 

\ 

exterior  ballistician. 

There  are  two  difficulties  in  extending  the  usual  methods  of  non¬ 
linear  mechanics  to  the  problem  of  missile  yawing  motion.  The  first 
difficulty  is  primarily  algebraic  and  the  second  conceptual.  Since  the 
yawing  motion  has  two  degrees  of  freedom,  the  problem  required  the 
solution  of  a  fourth  order  system  of  equations  and,  hence,  the  appli¬ 
cation  of  the  Kryloff-Bogolinboff5  techniques  can  require  rather  compli¬ 
cated  algebraic  operations.  Secondly,  vibrations  of  this  fourth  order 
system  can  take  on  two  frequencies  so  that  the  averaging  step  of  the 
K  -  B  method  would  require  an  averaging  of  two  different  frequencies. 
Although  this  is  mathematically  possible,  the  precise  meaning  of  such 
an  average  is  vague. 

In  this  report  it  will  be  shown  that  for  an  important  subset  of 
the  set  of  fourth  order  differential  systems,  the  fourth  order  equation 
in  a  real  dependent  variable  can  be  replaced  by  a  second  order  analytic 
equation  in  a  complex  variable.  This  use-  of  the  complex  variable 
introduces  important  simplifications  into  the  non-linear  problem.  In 
Reference  6  it  is  shown  that  with  the  proper  selection  of  geometrical 


The  angle  between  the  missile's  axis  and  the  tangent  to  its  trajectory 
is  called  the  yaw  angle.  The  yawing  motion  is  usually  described  by 
two  variables .  Two  common  choices  are  the  Eulerian  angles  locating 
the  missile  axis  in  wind-fixed  coordinates  and  the  direction  cosifles 
of  the  velocity  vector  with  respect  to  missile-fixed  coordinates. 
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variables , the  system  of  differential  equations  describing  the  yawing 
motion  can  be  a  member  of  this  subset.  In  this  report  it  will  be  shown 
that  for  a  large  class  of  non-linearities  the  averaging  step  of  the 
K-B  method  need  only  be  done  for  a  single  frequency  and  that  most  im¬ 
portant  aerodynamic  non-linearities  lie  in  this  class. 

By  means  of  an  "amplitude  plane"  the  results  of  the  analysis  are 
presented  in  a  simple  and  revealing  form.  It  is  shown  that  the  whole 
character  of  the  motion  can  be  described  when  the  location  and  nature  of 
a  small  number  of  singularities  are  determined.  The  effect  of  a  cubic 
Magnus  moment  is  discussed  in  some  detail.  Finally, predictions  for 
limit  motion  of  two  different  cases  of  quintic  Magnus  moments  are 
compared  with  exact  numerical  integrations  made  at  the  Naval  Proving 
Ground.^ 

Although  only  the  homogeneous  part  of  the  complex  yaw  equation 
is  treated  in  the  report  proper,  the  effect  of  gravity- induced  "yaw  of 
repose  and  the  effect  of  small  asymmetries  is  considered  in  an  appendix. 
The  direct  substitution  method  of  Reference  6  is  used  and  the  results 
are  compared  with  additional  Naval  Proving  Ground  calculations.® 

2 .  EPICYCLIC  SUBSET 

The  most  general  linear  homogeneous  fourth  order  system  of 
differential  equations  with  constant  coefficients  may  be  written  in 
the  form 

x"  =  c±  x'  +  c2  x  +  c3  y>  +  ck  y  (la) 

y"  =  c5  x»  +  c6  x  +  c?  y»  +  eg  y  (lb) 

where  c^  are  constants  and  primes  indicate  derivatives  with  respect  to 
the  independent  variable,  p.  If  Eq.  (lb)  is  multiplied  by  i  and  added  to 
Eq.  (la)  and  the  complex  variable  5  =  x  +  iy  is  introduced,  the  following 
complex  equation  may  be  written 

5"  +  A  V  +  B  |  +  C  V  +  D  1  =  0  (2) 

where  the  complex  coefficients,  A,  B,  C,  D,  are  linear  combinations  of  the 
ci •  W*1*211  c  and  D  vanish,  Eq.  (2)  becomes  an  analytic  differential  equation 
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In  the  complex  variable,  and  can  be  easily  solved.  Because  of  this 
simplicity  It  Is  of  some  Interest  to  determine  the  conditions  under 
which  it  Is  possible  to  transform  Eq.  (2)  to  this  special  form  by  means 
of  a  reversible  linear  trans format ion.  For  a  reason  which  we  will  give 
later  in  this  section  this  subset  of  the  set  of  all  differential  systems 
of  the  fourth  order  with  constant  coefficients  will  be  called  the 
epicyclic  subset. 

The  most  general  linear  transformation  can  be  written  in  the  form 

X  =  m^l  +  DVjI  (3) 

where  m^  and  m^  are  complex  constants.  The  inverse  equation  can  be 
'obtained  by  eliminating  T  between  Eq.  (3)  and  tba.jagajugafe  ■  of  (3). 


I  = 


X 


“2 


w 


nLjm  - 

Thus  the  restriction  to  reversible  transformations  is  equivalent  to  the 
relation 

ml^l  “  ^  0  (5) 

Without  loss  of  generality  we  can  assume  that  m^  is  non- zero*,  and,  hence, 
relations  (3)  and  (5)  become 


i© 


X  =  me  (|  +  si) 


m 


(1  -  ss)  £  0 


±0 


where  me  =  m. 


V 


(6a) 

((6b) 


s  =  — ,  and 
ml 


m  is  real. 
-10. 


I 


10  — 
x  -  se  x 


m(l  -  ss) 


(7) 


*  . 
If  m^  were  zero, 


we  would  consider  the  conjugates  of  Eq. 


(2)  and  (3). 
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Substituting  Eq.  (7)  in  Eq.  (2), 

X"  +  (A  -  C  s)x'  +  (3  -  D  I)x 


2i6 


[sr- 


'+  (As  -  c)X’  +  (Bs  -  D)X  I  =  0 

If  X,"  is  eliminated,  between  Eq.  (8)  and  its  conjugate, 

(1  -  ss)x"  +  (A  -  Cs  +  Cs  -  A  ss)xf 


+  (B  -  Ds  +  Ds  -  B  ss)x  -  e 

+  [5  s2  +  (B  -  B)s  -  D  IJ  =  0 
where  (1  -  ss)  is  non  zero  by  equation  (6b). 


! 


2i©  JU  2 


C  s  +  (A  -  A)s 


(6) 


(9) 


The  requirement  that  Eq.  (2)  belong  to  the  epicyclic  subset  means 
that  it  must  be  possible  to  select  s  so  that  the  coefficients  of  X’  and  X 
vanish.  From  Eq.  (9)  it  can  be  seen  that  m  and  ©  cannot  make  the  coef¬ 
ficients  of  X’  and  X  vanish  and  so  they  may  be  arbitrarily  selected.  If 
we  make  the  definitions  A  =  A1  +  iAg,  B  =  B1  +  iB2,  C  =  |c|eior,  and 

D  =  I  D  I  e^°,  then  we  have  the  following  pair  of  quadratic  equations  for  s. 


c 

s2  +  2iA2seitX  - 

C 

D 

s2  +  2iB2sek°  - 

D 

e2ia  =  0 


2io)  _ 
e  =  0 


(10) 

(ID 


If  neither  C  no.rD  vanish,  the  solutions  of  Eqs.  (10  -  11 )  are: 


A2  1  / 

-2 - 

s  * 

c 

-  1 

n  J 

2 

2  J 

B2 

PI  ^ 

D 

“  1 

Ka  -  |) 


K®  -  5) 


(12) 


(13) 


Since  s  is  unity  when  either 


A2 

B2 

c  or 

D 

is  less  than  or  equal  to  one 


and  this  is  contrary  to  relation  (6b),  we  have  the  restrictions: 

71  (14) 


*2 

B2 

C 

71, 

D 

10 


I 


Under  restriction  (l4), srjcan  satisfy  both  Eq.  (12)  and  Eq.  (15)  only  if 


B, 


_^2  ___ 
i°r  idi 

in  the  single  equation 


2  *2 
and  q  =  cjj  or  .-g-  = 


-B_ 


and  a  =  o>  +  rt.  This  can  be  given 


^2 

C 


B„ 


(15) 


By  means  of  Eqs.  (10  -  15)  we  can  now  state  the  following  theorem: 

Theorem  Eq.  (2)  is  a  member  of  the  epicyclic  subset  if  one 
of  the  following  conditions  is  satisfied: 

(1)  C  =  D  =  0 

I  I 

7  1 


(2)  0  =  0,^  =  0, 


(5)  D  =  0,  B2  =  0, 


C 


>1 


Ag  B, 
C 


(*>  4-4 


2 

D  ’ 


c 


7  1. 


For  the  trivial  case  (l),  s  is  zero;  for  the  other  cases  it  is  fixed 
by  either  Eq.  (12)  or  Eq.  (15)* 

By  means  of  such  a  linear  transformation  all  equations  belonging 
to  the  epicyclic  subset  may  be  written  in  the  form 


+  (ax  +  i&2)Xf  +  (^  +  ib2)\  =  0 


(16) 

(t  a  +  10*  )p 

where  a^,  a2,  b^,  and  b2  are  constants.  Substituting  \  =  e'' 

in  Eq.  (l6)  and  separating  real  and  imaginary  parts,  we  see  that  0*  and 
a  must  satisfy  the  following  equations: 


fj.2  ml  p  •  i,  .  a(a^  -  -  3 


a  = 


a  0*  +  t 
2  0*  +  a. 


(18) 
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Under  the  usual  assumption  of  small  damping  during  a  cycle,  the  a  term 
in  Eq.  (17 )  can  be  neglected  and  it  reduces  to 


0* 


d 


2 


J 


1,2. 


(19) 


+  4  b1>  0 


If  the  above  inequality  is  not  satisfied  or  the  damping  is  not  i  over 
a  cycle,  Eqs.  (17  -  18)  lead  to  the  need  for  the  solution  of  a  fourth  order 
equation.  This  difficulty  can  be  avoided  by  not  separating  Eq.  (16)  into 
real  and  imaginary  parts  after  making  the  substitution. 


All  the  cases  treated  by  the  methods  of  this  report  will  be  restricted 
to  those  with  small  damping  over  a  cycle  and  Eq,  (19)  will  apply. 

Since  a 2  =  -  (01*  +  02')>  symmetric  forms  of  the  equation  for  the 
damping  exponents  can  be  obtained  from  Eq.  (18). 


a. 


h  '  +  *2 
K'  -  A,' 


(20) 


au 


al  02*  +  ^ 

02*  -  0-l* 


(21) 


The  equation  for  the  general  solution  to  Eq.  (l6),  therefore,  is 


\  =  K, 


i(01O  +  0!*P) 


1(0. 


20 


02f 


P) 


(22) 


-a,p 

where  Kj  =  KJQ  e  ;  0  j '  and  0^  are  given  by  Eqs.  (19-21);  and  0JQ 

and  KjQ  are  constants.  This  solution  is  a  linear  combination  of  two 

complex  vectors  which  are  rotating  at  certain  fixed  frequencies  and  are 
exponentially -damped.  Since  the  curve  swept  out  in  the  complex  plane  is 
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called  an  epicycle,  the  reason  for  the  name  of  this  subset  is  now  clear. 

The  solution  for  i  when  Eq.  (2)  belongs  to  the  epicyclic  subset  can 
be  obtained  by  use  of  Eq.  (4).  s  is  determined  to  be  a  solution  of  Eqs.  (10 
11 ),  m  and  0  are  arbitrary  and  for  convenience  m  will  be  made  unity  and  0 
will  be  fixed  at  zero. 


1  -  ss 

i0-i  i02  -101  'i02 

e  +  K2  e  +sKxe  tsy 


1  -  ss 


(23) 


where  0^  =  0jO  +  0  Jp. 

3.  SOLUTION  OF  THE  NON-LINEAR  EPICYCLIC  EQUATION 

Although  the  treatment  of  non-linear  fourth  order  systems  is  usually 

quite  laborious,  the  procedure  for  systems  which  may  be  linearized  to 

members  of  the  epicyclic  subset  is  much  simpler.  In  this  section  the 

approximate  solution  to  equations  of  the  form  of  Eq.  (l6)  but  with 

coefficients  a  ,  b^  b2,  which  are  functions  of  X,*,  and  \* , 

will  be  considered.  The  method  of  solution  which  will  be  employed  will 

5 

be  essentially  that  of  Kryloff  and  Bogdilljibfbff , 


This  method  is  based  on  a  pertubation  of  the  solution  of  the 
linear  equation  with  no  damping.  If  this  solution  is  identified  by  X, 
then 


where 


~  =  K10  e 


10O 


*20 


(24) 


0j  -  0j o  -  0jP 


rJ 

^  *a2 
*r  — 


fa  ,  ~~p 
+  j  a2  +  4  b1 


^  +  4^  7  0 

a  ,  b1  constant  parts  of  a£  and  b±  (their  values  when 

X.*  =  X.  =  o. );  and 

0.~  are  constants. 

JO'  rj0 
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The  complete  non-linear  differential  equation  can  be  written  in  the  form: 


\  +  is.^  \ f  +  b^  \  — 


[al  +  ^a2  -  a2}J  "  |^bl  '  V  + 


=  f(X,  X,  X',  X1). 

where  a_L,  ag,  b^  bg  are  functions  of  X,  X,  X1 ,  X’. 


(25) 


The  pertubation  process  is  quite  similar  to  the  well-known  method 

of  variation  of  parameters.  First  the  solution  is  assumed  to  be  that 

expressed  by  Eq.  (24)  with  the  constants  Kjq,  replaced  by  unknown 

functions  of  p,  K.(p)  and  \|r  (p). 

i(^  +  0j)  i(+2  +  \) 

+  1^6 


.  •  X  =  e 


(26) 


This  equation  is  then  differentiated  to  yield: 
X1 


y  i(t  +  0  )  y  i (*„  +  %) 

.  10^  Kx  e  1  1  +  102'  e  2  2 


I 


,  f  !(*  +?)  i(^p  +  0P? 

+  (K^  +  +  (Kg*  +  i*2'K2)e  2  2  (27) 

The  last  two  expressions  in  Eq.  (27)  are  set  equal  to  zero  so  that 


,  .  i(^  +?)  i(i|r  ?+0p) 

(K-l*  +  il|f1,K1)e  +  (Kg*  +  i-VK2)e  2 

and  Eq.  (27)  then  simplifies  to 


=  0  (28) 


• 

X'  =  10  ^  e 


i(1J,1  +  0-|_)  ~  _  i(t2  +  02) 


+  i  02*  K2  e 


Differentiating  again  we  see  that 

,2  i(t1  +  0X)  y  g  i(*  +  0  ) 

1  K,  e  -  0  '  Kr  «  ■  2  2 


X"  =  -  0X*  Kx  e 


iU]  +  K)  ~ 

)e  +  i  0o’(K»  +  i*  'I 


+  1  0i,(K1*  +  i^Kx)e  x  x  +  i  02»(K2*  +  iHfg'Kg) 

If  Eqs.  (26,  29,  30)  are  substituted  in  Eq.  (25), the  terms  without 


(29) 


(30). 


i(*2+02) 
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derivatives  of  the  parameters  must  satisfy  the  homogeneous  equation  and 
therefore  cancel. 


^’(V 


i(VHO  ^ 

i+jK^e  1  1  +  i02,(K2*  +  i+2,K2)e  =  f(p) 


(31) 


(Kg*  +  it2 *K2 )  can  he  eliminated  between  Eqs.  (28)  and  (3l). 


-i(+!  +  0X) 


rj 

=  -  (0 


i(^x*  -  ^)K± 

t  ill  ,  ^  . 

,i'  .  02»)  "  J  -  &i  +  i(a2  -  a2) 


(32) 


& 


[b2  -  -  b1) 

.  _  rv> 

where  0  =  02  +  ^  -  0X  - 


i +  %  f 


v  ♦  k  ^  *10 


This  simple  derivation  of  Eq.  (32)  shows  the  considerable  reduction  in 
algebraic  complication  which  the  use  of  the  complex  variable  and  the 
restriction  to  the  eptcyclic  subset  have  introduced. 


At  this  point  in  the  K-B  method,  due  to  the  difficulty  in  solving 

the  differential  equations  for  the  parametric  functions,  averages  of  the 

right  sides  of  the  parametric  equations  are  taken.  In  general,  the 

right  hand  side  of  Eq,  (32)  has  two  basic  frequencies  and  the  necessary 

averages  have  a  rather  vague  physical  meaning.  If  only  the  difference 

frequency  0)  is  allowed  in  the  functional  dependence  of  a^,  a2,  b^,  and 

b2  on  X,  X,  X.',  and  X',  a  single  average  would  be  needed  and  this  average 

would  have  a  simple  meaning.  This  means  that  only  those  non-linearities 

15 

which  can  be  associated  with  a  body  of  revolution  are  considered.  For 
example,  under  this  restriction  there  are  four  possible  quadratic  com¬ 
binations  of  X,  X,  X*  and  X1:* 


*  The  effect  of  damping  has  been  neglected  in  the  calculation  of  X*  and  X*  . 
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62  =  U  =  K12  +  K/  +  2K  X_  cob  (0) 

X'X*  =  +  02*2  K2  +  2  01»02*  KxK2  cos  {0, 

VK12  +  4,K22  +  KiK2^l*  +  02*  e_i^) 

0l'Ki2  +  02»K22  +  K1K2(01*  e'^  +  02*  e^) 


U,e  -  .1 


X'X  -  i 


(33) 

(34) 

(35) 

(36) 


From  Eq.  (35  -  34)  we  Bee  that  both  the  magnitude  ...of.  X  and  the  magnitude  of 
its  derivative  have  the  frequency  of  .0 .  and,.  hence,  the  averaging,  .will  be 
made  over  a  period  of  these  amplitudes . 

It  is  now  necessary  to  assume  that  the  modal  amplitudes,  Kj,  change 
slowly  over  a  period  of  the  difference  frequency.  With  this  assumption 
in  mind  we  can  average  the  right  side  of  Eq.  (32)  over  a  period  of  0  and 
separate  the  result  into  real  and  imaginary  parts ^ 


V 

K 


-1 


i  2*(0'1*-5s'2') 


2  it 

f,f2 
2  K1 

r 

- 

a-L  J^1  +  5 

cos  0 

+  h2 

IT  k2 

t  ^ 

cos  0 

^0 

(ap  -  a2)$2»  -  (b 


~*1  K? 

>i-bij  icf 


■  in  0 


[  d0  =  -  a±  (K.^)  (37) 


2jf 


*1*  = 


2m  (Jf. 


i-4’)  r 


(a2  -  a2) 


[W2*| 


cos 


0 


-  (b,  -  bx) 


r  k2 

rsS 

1  +  rp-  COS  0 
ft.-. 

- 

al02*  +  b2 

L  l  _ 

—  _ 

K  1 

K"  sin  0  > d0  .  (38) 
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From  symmetry,  the  equations  for  the  second  mode  of  oscillation  are 


Thus  for  each  cycle  of  0,  the  K-B  approximation  is  a  damped  epicycle  with 

damping  exponents,  a.,  given  by  Eqs.  (37)  and  (39)  and  with  frequencies, 

J 

0  *  equal  to  0'  +  where  the  i|r  *  *s  are  given  by  Eqs*  (38)  and  (40). 
J  J  J  J 

The  damping  exponents  predicted  by  Eqs.  (37)  and  (39)  for  the  linear  case 
are  precisely  those  given  by  Eqs.  (20  -  21 ). 

If  a2  and  b.^  had  been  chosen  to  be  the  average  values  of  a2  and  b^ 
instead  of  their  values  for  zero  amplitude  motion,  the  expressions  for  \|r  1 
would  have  been  simplified.  The  implications  of  such  a  choice  are 
described  in  Appendix  C.  In  most  cases,  however,  the  advantages  of  this 
choice  do  not  outweigh  the  alfeebraic  complexities  which  are  introduced. 
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.  4.  APPLICATION  TO  MISSILE  YAWING  MOTION 

With  the  aid  of  the  ipreaedlng;  two  sections  on  the  approximate 
solution  of  non-linear  equations  of  the  epi cyclic  type,  the  problem  of 
honkliaea&x1  yawing  motion  of  a  symmetric  missile  in  free  flight  can  be 
easily  handled.  For  this  case  the  complex  plane  is  the  planelperpen- 
dicular  to  the  missile’s  axis.  The  complex  vector,  lies  on  the 
intersection  of  this  plane  with  the  plane  of  the  yaw  angle  with  magni¬ 
tude  equal  to  the  sine  of  that  angle.  In  Reference  6  it  is  shown  that 
the  general  equation  of  yawing  motion  in  a  non-rotating  coordinate 
system  is 


X"  + 


H  +  Jg  -  -  iv  XtJ  + 

-  M  +  ij  1 

J_r 

-  iv  T 

-4 

_ 

where 


H  m  j^L  ^  +  k2  (KH  "  ^SlA^ 

!>  =/l  -  &2  (cosine  of  yaw  angle) 

&  =  (magnitude  of  sine  of  yaw  angle) 


\  =  G  - 


£' 

.IT 


v  = 


CD-^d 


u 


v  -££ 

L  m 


'djr  (62) 

X  -  (i2r 


L_d6 


G  =  f '  - 


m  ^  "  k2  KH^  ‘  Jg  +  iv 


(41) 
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(go  +  igJd 

r  =  — - ^ -  -  J  X 

h  2  g 


u 

axial  radius  of  gyration  in  calibers 
kg  transverse  radius  of  gyration  in  calibers 
p  density  of  air 
d  diameter 


m  mass 


o>1  axial  spin 
u  magnitude  of  velocity 


g^p  trajectory  component  of  the  gravitational  acceleration 

gg,  gj  components  of  gravitational  acceleration  perpendicular  to 
missile's  axis  and  the  K^'s  are  aerodynamic  coefficients  which  are 
defined  in  the  Table  of  Symbols.  Since  only  the  homogeneous  equation 
has  been  considered  in  Sections  2  and  3>  "the  discussion  of  the  inhomo¬ 
geneous  part  of  Eq.  (41)  will  be  deferred  to  Appendix  A. 


The  most  common  non-linear  forms  of  Eq.  (4l)  arise  from  a 

2 

dependence  of  the  aerodynamic  coefficients  on  5  .  Since  any  function 
2 

of  5  is  an  even  periodic  function  of  0^  -  02*  its  Fourier  series 
expansion  is  a  cosine  series  and,  hence,  most  of  the  sine  terms  in 
Eqs.  (37  ~  40)  vanish.*  Eqs.  (37  -  40)  can  be  written  as 


+1 


i 


*2 


i 


-1 


2« 


2«(?1'-^2')  J 


-1 

2«(?2»  -  ?!*) 


Kp 

—sin  0 
K1 


r  d0 


*  T 


(43) 


01*  j^sin  0 


* 

Since  lx  and  J  '  are  the  derivatives  of  even  functions,  they  are  odd 

J_j 

functions  and  so  have  sine  expansions.  Fbr  theses  quantities  oiilyi^the 
s  ine^,  terms.;  haye,  non-hero  contributions. ,  ...  •.  ’  ■  ’.  ... 
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V 


-1 


2rt 


K, 


2  2*(?2«-^') 


H(62)+Je 


-  ^  K  ”1  ""  K 

^'^i'  ir  cos  _  ^  1  +  ircos  0 


K 


+  £JL  Kj  sin  0  ^  d0  =  -  0^,  Kg) 


(^5) 


where 


V  -s  [»i  / 


v2  -  4m  7  0 


v2  -  4m 


and 


Since 


2* 


51  J  l1  8ln  *  **  -  S 


2  it 

r 


t1  -  V '  sin  0  d0 
2(1  -  52) 


2« 


55  KxK2  sin20  (0^  -  02*)(1  +  S2 )  d0 


(01,-?2')  KiK2  f1  +  Ki2  +  */) 


y  and 


*  —2 

The  inequality  v  -  4m  y  0  is  equivalent  to  the  requirement  that  the  missile 
be  gyroscopically  stable. 9  For  positive  spin,  0  *  is  the  larger  frequency 
and  is  usually  called  the  nutational  frequency  while  0’,  the  smaller 
frequency,  is  called  the  precessional  frequency.  On  page  19  of  Ref.  9, 
it  is  shown  that  this  use  of  the  terms  "precession"  and  ^nutation"  is  not 
compatible  with  their  use  in  the  theory  of  the  top. 
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-v +  V 


then 


0!1  ;  ^  ~ 


a,  02f  4  <i+4  +  4) 


2« 


2*(£L‘-?2'); 


■J 

M(52)-M(0) 

"0 

_  — 

^  4  (i+4+k2j 

0  •  =  ^  -  -i-± - + 


2rt 


2rt(0j_  -  fg) 


M(B  )-M(0) 


0 


(47) 

K2 

(l  +  — cos0)d.0 
K1 

(48) 

K1 

(1  +  pr-cos0)d0 


Although  Eqs .  (47  -  48)  for  the  equivalent  linear  frequencies  are  of 
some  interest,  the  primary  concern  of  a  designer  is  the  behavior  of  the 
amplitude  of  the  motion.  For  the  type  of  non-linearities  which  are  usually 
encountered,  it  will  be  convenient  tb  ddecribife.-tMi  motion,  lih  ,tferms  of. 
the  squared  amplitudes  of  each  mode.  In  the  more  detailed  analysis  that 
follows  we  will  make  the  convenient  but  not  necessary  assumptions  of 
linear  damping  moment  and  linear  lift  force  (H  =  constant,  J^’  =  0). 

Eqs .  (44  -45)  then  reduce  to 


(K-,2)'  ?  2 

— ±5 -  =  -  2a1  (Kx  ,  Kg  )  (49) 

K1 


(K2)'  2  2, 

- -2a2  (Kx  ,  Kg  ) 

K2 


<H  +  Jg$i'  -  s  1 

2n 

^  2 
l  v  T(6  ) 

0 

K2 

1+77-  COS  0 
_  1 

d  0 

<H  +  V4'  •  h  \ 

V  -  02' 

i  v  T(52) 

'  0 

”  K1 

1  +  rf-  COS  0 

L  K2  J 

d0 

02  ■*  -  iV 

(50) 


(51) 


(52) 


These  frequencies  are  measured  in  the  non-rotating  coordinate  system. 
The  relations  between  the  non- rotating  coordinate  system  and  the  fixed 
plane  coordinates  sire  given  in  Reference  6. 

In  Appendix  D  the  interesting  special  case  of  a  non-linear  damping 
moment  and  zero  spin  is  considered. 
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where 


^  =  K^2  +  K22  +  008 

t 

5.  THE  AMPLITUDE  PLANE 

In  Reference  6  the  analysis  of  non-linear  spark  range  data  was 
developed  in  some  detail.  The  yawing  motion  over  the  relatively  short 
flat  trajectory  under  observation  was  assumed  to  be  epicyclic  with 
constant  frequencies  and  damping  exponents  which  were  related  to  the 
appropriate  effective  magnitudes  of  yaw.  For  the  cubic  Mggnus  and 
static  moments,  which  were  assumed,  Eqs .  (47  -  48)  and  (51  -  52) 
reduced  to* 


*The  quart ic  terms  in  Eqs.  <47  -  48)  are  omitted  for  B  <£  sin  15°. 
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where 


n 


M 


-  Z>  k2'2  ^  b21C  4  “o  +  ^ 


k=0 


+  B 


n 


-.z 


i  J  -  k  '  J 

ij2k  ^  ^Pk 

k=0  L  - 


pv  2 

6  =  TQ  +  T2  5  . 


.  JJ  =  —  K.  . 

1  m  a 


Using  Eqs.  (53  -  56)  a  large  number  of  spark  range  firings  were  analysed 
with  outstanding  success.  The  swerving  motion  was  treated  in  a  similar 
fashion  and  cubic  lift  force  and  cubic  Magnus  force  coefficients  were 
obtained  which  showed  excellent  consistency  with  their  corresponding 
moment  coefficients. 

Although  this  work  on  spark  range  firings  was  so  successful,  it 

suffered  from  its  limitations  to  cubic  non-linearities  and  short  portions 

of  trajectories.  In  Reference  6,  the  influence  of  polynomiaLlexpemsion 

of  M  and  T  was  calculated  up  to  l4th  degree  polynomials*  and  Eqs.  (47  - 

2 

48,  51  -  52)  will  apply  when  M  and  T  are  arbitrary  functions**  of  5  . 

Thus  the  basic  limitation  is  the  restriction  to  short  trajectories. 


*  The  effective  yaws  of  that  report  are  related  to  the  integrals  of 
this  report  by  the  definitions 


In  Appendix  B  effective  values  of  negative  powers  of  B  are  computed 
and  a  possible  use  indicated. 
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This  restriction  lies  In  the  assumption  that  the  yawing  motion  is 
epicyclic  with  constant  frequencies  and  constant  damping  exponents. 
(These  constants  depend  on  the  effective  amplitude  of  the  motion)} 

A  study  of  Eqs.  (49  -  50)  shows  how  this  restriction  can  be  relaxed. 
The  theory  requires  that  these  quantities  be  constant  over  a  period  of 
6  but  they  may  vary  over  longer  intervals.  Dividing  Eq.  (50)  by  Eq. 

9 ),  we  can  write  a  single  first  order  non-linear  equation  for  this 
variation 


(K22).»  d(K/)  K. 22  a2  (Kj2,  k/) 

(K^)'  (K2,  K^) 


(57) 


Equation  57  describes  the  character  of  the  yawing  motion  by  means 
of  the  movement  of  a  point  in  the  K±  ,  plane  which  we  will  call  the 
amplitude  plane.  For  any  point  in  this  plane  the  equivalent  linear 
frequencies  and  damping  exponents  can  be  calculated  and  thus,  except  for 
the  phase  angles,  0jq,  the  motion  is  completely  determined.  As  will  be 
seen,  this  amplitude  plane  will  have  a  number  of  similarities  with  the 
phase  plane  associated  with  the  one  degree  of  freedom  problem.  Although 
the  usual  four  dimensional  phase  space  associated  with  two 'degrees  of 
freedom;  reduces  to  three  essential  dimensions  for  the  epicyclic  subset, 
it  will  not  in  general  reduce  to  two  .dimensions  and,  hence,  the  amplitude 
plane  is  definitely  not  a  phase  plane. 

Differential  equations  of  the  form  of  Eq.  (57)  have  been  treated  in 
some  detail  by  Poincare10  who  showed  that  the  essential  properties  of 
the  solution  curves  are  fixed  by  the  location  and  type  of ffrrir  singularities. 
Singularities  are  points  for  which  both  the  numerator  and  the  denominator 
of  the  right  side  of  Eq.  (57)  vanish.  They  can  vanish  in  four  different 
ways..  These  ways  correspond  to  the  following  points! 

(1)  the  origin; 

2 

(2)  the  intercepts  of  the  =  0  curve; 

L  1  2 

(3)  the  K2  intercepts  of  the  =0  curve;  and 

(4)  the  intersections  of  the  zero  damping  curves  =  0 

and  a2  =  0. 
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According  to  Poincare's  classification  there  are  four  possible  types 
of  first  order  singularities:  nodes,  saddles,  spirals,  and  centers.  In 
Fig.  1, a  node  and  a  saddle  are  shown,  while  a  spiral  appears  in  Fig.  3b 
and  a  center  in  Fig.  4.  The  particular  Curve  along  which  the  actual 
yawing  motion  moves  is  determined  by  the  initial  modal  amplitudes  as 
specified  by  the  initial  conditions.  If  the  coordinates  are  translated 
so  that  the  singularity  is  at  the  origin,  Eq.  (57)  has  the  form: 


dy  _  ax  +  by  +  P(x,  y)  (=3 

dx  cx  +  dy  +  Q(x,  y) 

where  ad  -  be  ^  0  and  P  and  Q  vanish  to  at  least  the  second  order  at 

the  origin.  The  criteria!  JFor  the  type  of  the  singularity  can  now  be 
stated  in  terms  of  the  coefficients  a,  b,  c,  d  and  their  discriminant 
D  =  (b  -  c)2  +  4  ad.  (See  page  kk  of  Ref.  11.) 

I.  The  singularity  is  a  node  if  (l)  D  ?  0  and  ad  -  be  <.  0  or 

(2)  D  =  0.  (Note  that  if  ad  =  0,  it  is  a  node  if  bh > 

II.  The  singularity  is  a  saddle  if  D  7  0  and  ad  -  be  7  0. 

III.  The  singularity  is  a  spiral  if  D  ^  0  and  b  +  c  ^  0. 

IV.  The  singularity  may  be  a  center  if  D  40  and  b  +  c  =  0; 

otherwise  it  is  a  spiral.  The  higher  order  terms  P  and  Q  must  be 

Ik 

considered  for  a  final  determination. 


For  Eq.  (57)>  the  character  of  the  singularity  at  the  origin  can  be 
i&Sily  determined  since  a  and  d  both  vanish  and  b  and  c  are  the  linear 
values  of  the  damping  exponents  a2Q,  o^q.  The  origin  is  a  node  when 
these  damping  exponents  are  of  the  same  sign  (a-j_oa20  7  ^  6X1(1  a  sati<ile 
when  they  are  of  opposite  sign  (o^o^O  ^  0)*  In  Fig.  la,  the  first 


possibility  is  shown  for  a1Q 


a 


10 


X 


a, 


20 


and  in  Fig.  lb,  the  second  for 


a, 


20 


is  shown.  The  direction  of  motion  is  not  shown  since  this 
depends  on  the  actual  sign  of  the  damping.  If,  for  example,  both  exponents 
are  positive,  the  origin  is  a  stable  node  and  small  amplitude  motion  must 
damp  to  zero. 


* 

If  ad  -  be  =  0  and  the  origin  is  a  singularity,  then  the  singularity  is 
at  least  second  order  and  the  P  and  Q  functions  must  be  considered. 
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2  2 

The  singularities  on  the  or  axes  can  be  located  by  consideration 
of  pure  mode  motion.  They  are,  however,  very  important  for  "almost"  pure 
modes  as  well.  For  pure  nutational  motion,  i.e.,  (l^2  =  0),  Eq.  (51)  reduces 
to 


^(k/,  0)  = 


(H+Jg)0^*  -  V  T(KX2) 
-  %' 


(59) 


Similarly  for  pure  precessional  motion. 


a9(o,  V)  = 


2  (H+J  )02*  -  v  TCk/) 


(60) 


Eqs.  (59  -  60)  could  have  been  obtained  by  taking  the  expressions  for 
linear  damping  and  replacing  the  constant  Magnus  moment  coefficient  by 
its  actual  non-linear  dependence  on  &2.  This  is  the  process  used  in 
Refs.  1  and  2  in  their  treatment  of  pure  mode  motions. 


In  order  to  consider  motion  near  a  pure  mode  of  amplitude  K, 
)  must  be  expanded  in  Taylor  series  ’about  this  amplitude . 


T(62)  =  TtK2)  +  T* (K2 )(B2  -  K2) 


(61) 


Substituting  Eq.  (6l)  in  Eqs.  (51  -  52), 


Ko2)  = 


(H+Jg)^!  -  v  TOC2)  -  vT*  (K2)  [Kl2  f  21C,2  -  K2^ 

V  ■ 


(62) 


2  2,  (H+Jg)iV-  -  vT*  (K2)  fe2  +  2KX2  -  K2] 

\T\' 


^(Kj.  ,Kg  ) 


(63) 


where  the  point  (K^2,  )  is  near  either  (K2,  0)  or  (0,  K2).  If 

singularities  on  the  nutational  axis  ( K ^  =  0)  are  considered,  their  location 

iayr.becDbtained  bly L s ettlng/ Eqo t  ( 59 j  :dqual).5^  eg roqr o  I dehiiif JitigL isachu aui  a 
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singular  point  by  the  coordinates  (YT ,  0),  Eq.  (59)  becomes 


(H+J -  v  T(K‘') 


-3.  =  0 


(64) 


If  Eq.  (64)  is  substituted  In  Eqs.  (62  -  65), 
(k/,  k/)  =  -  A  [(Kl2  -  K2)  +  2Kg2 


CLr 


(I if,  Kg2)  =  |H,+Jg  +  AK2  +  A  ^Kg2  +  2(KX2  -  K2) 


(65) 


(^6) 


where  A  =  .V  •  Substituting  Eqs.  (65  -  66)  in  Eq.  (57)  and.  shifting 

“  V 

2  2  2 

the  origin  to  the  singular  point  by  the  translation  x  =  -  It,  y  =  Kg  , 

we  see  that 


and 


dy  ,y  liH+Jg+A^) 

+  A(y  +  2x)J 

^  (x+K2)  [- 

A(x  +  2y)J 

. ’.  a  =  0 

c  *-AK' 

b  =  H+J  •+  AK2 

d  =  -  : 

g 

D  =  +-2AK2)2 

g 

^0 

2AK2, 


(67i) 

(68) 

(69) 


According  to  relation  (69)  the  singularity  must  be  either  a  node  or  a 
saddle.  It  will  be  a  node  if  A(H+Jg  +  AK2 )  is  negative  and  a  saddle 
otherwise.  Conversely,  if  the  singularity  Is  on  the  precessional  axis, 
it  can  be  located  by  setting  Eq.  (60)  equal  to  zero.  Since  Eqs.  (59  -60 ) 
are  symmetric  in  the  subscripts  1  and  2,  the  test  for  type  of  singularity 
can  be  obtained  by  Interchanging  subscripts  in  the  above  discussion. 
Therefore,  the  precessional  singularity  will  be  a  node  if  -A  (H+Jg  H-  AK2 ) 
Is  negative  and  a  saddle  otherwise. 
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The  situation  for  singularities  which  are  not  on  the  coordinate  axes 
is  much  more  complex.  In  order  to  treat  these  singularities  it  is  neces¬ 
sary  to  know  much  more  than  the  value  of  T  and  T'  at  a  point.  To  illus¬ 
trate  the  proper  procedure  for  handling  these  singularities  the  possible 
amplitude  planes  for  a  cubic  Magnus  moment  will  be  discussed  in  some  detail. 

According  to  Eqs.  (55  -  %)  the  curves  of  zero  damping  are  lined  des¬ 
cribed  by  the  following  equations: 


a. 


10 


+  a!2  (K1  +  ^2  >  =  0 


(70) 


a20  "  A12  (K2  +  ^1  )  "  9 


(71) 


where  a. 


K 


10 


H  0  »  -  v  T, 


-  % 


H  0, 


V  T, 


a, 


'20 


f£>*  -  \ 


0 


ai2  _  ~  a22 


-  V  Tr 


The  Intersection  of  the  lines  is  a  singularity  and  It  has  coordinates 


riO  +  2a20  a20  +  2aiON'| 

- t-t -  J.  In  addition  to  this  there  are  singularities 

'**12  / 


3a. 


12 


at  the  origin  and  the  points 


/  O10 


,  o)  ,  fo,  These  three 


\  ai2  /  V  '■*12/ 
singularities  cap,  however,  be  treated  by  methods  already  developed. 

Since  only  the  first  quadrant  of  thelamplltude upland  h&d-'phyilcd&n&eanlng, 
the  intersection  of  the  zero  damping  lines  has  Importance  only  when  it 
lies  in  the  first  quadrant. 


O10  +  2a20  a20  2ai0  _  „ 

z  ^  X  O 


3  a. 


12.. 


-  3  a. 


12 


(72) 
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9  a. 


If  these  are  multiplied  together  and  divided  by 


12 


,  then 


2  a, 


'20 


'a. 


+  21/5  +  -^\o 

a20  A  2  20  7 


10 


•  •  a  ^  2 

a20 


(73) 


(74) 


Comparing  inequality  (74)  with  the  original  inequalities  (72),  we  see  that 
the  intersection  will  be  in  the  first  quadrant  when  _ lies  between 


a, 


20 


-  2  and  -  ^  and  a2Q  has  the  same  sign  as  a^g. 


In  order  to  classify  this  singularity  Eqs.  (55  -  56)  are  now  placed 

.s  ti 
LO  + 

Tal 


in  Eq.  (57)  and  the  origin  is  translated  to  the  point  of  intersection  by 

2  Q10  +  2  a20  2  a20  +  2  a10 
the  translation  x  =  K-,  -  — 5— -  ,  y  =  Kg  + 


*12 


TS 


12 


Jjy  («20  +  2  a10^y  +  2x)  “  5  a19  (y  +  2x)y 
*  *  dx  = 


12 


(a1Q  +  2  a20)(x  +  2y)  +  5  a12  (x  +  2y)x 
The  coefficients  of  Eq.  (5®)  can  now  be  obtained  from  Eq.  (75)* 

C  "  Q10  +  2  °20 

i  .  2  (a10  +  2  a2Q) 


(75) 


E  -  2(“20  +  2  “lO5 


(76) 


b  =  +  2 

20  10 


.  .  D  =  53  a10  +  78  aioa20  +  ^  a20 
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a202 


10 

20 


+  1.812 


aio 

+  *552 

a20 


(77) 


d  -  be  =  3(a10  +  2  aon)(aon  +  2  ain) 


’20''  20 


‘10 ; 


b  +  c  =  3(a10  +  °20) 


(78) 

(79) 
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If  the  point  of  intersection  lies  in  the  first  quadrant,  be  -  ad  must  be 
always  positive.  Applying  the  criteria  for  singularities  to  relations 
(72,  75  -  11),  we  can  make  the  following  statements  for  the  singularity 
at  the  point  of  intersection. 

0*1  ,  QL 

(1)  It  is  a  node  if  -  2  <— —  =  -  1.812  or  -  .552^  —  .5. 

a20  a20 

(2)  It  is  a  spiral  if  -  1.812  <  -  .552  and  —  £  -  1. 

a20  a20 

?10 

(5)  It  may  be  a  center  if  ^ —  =  -  1.  (An  application  of  the  criteria 

u20 

of  Ref.  14  shows  that  it  is  a  center.) 


In  Figures  1  -  k*,  all  but  one  of  the  possible  different  amplitude 
planes  for  a  cubic  Magnus  moment  have  been  drawn  by  the  EBL  Analog 
Computer.  If  a2_oa20  or^8^n  Is  a  node  and  only  one  of  the  lines 

of  zero  damping  falls  in  the  first  quadrant.  Since  the  intercept  of  this 
line  has  to  be  a  saddle  point,  the  amplitude  plane  is  that  shown  in  Fig.  2. 
This  figure  shows  a  characteristic  property  of  non-linear  equations, 
namely  a  dependence  of  the  form  of  the  solution  on  initial  conditions.  If 
the  origin  is  a  stable  node,  the  motion  is  down  toward  the  -  axis  and, 
hence,  for  certain  initial  conditions,  the  motion  goes  to  zero  while  for 

p 

others  it  goes  to  large  values  of  K,  . 

If  Ctjo a20  ^  ttie  0riSln  is  a  saddle  and  either  both  lines  go 

through  the  first  quadrant  or  neither  doe#!  .'When  both  lines  go  through 

the  first  quadrant  but  do  not  intersect  there  (Fig.  5a),  our  criteria 

reveal  that  the  singularity  on  the  line  closer  to  the  origin  is  a  node 

while  the  other  singularity  is  a  daddle  point.  Thus,  for  a  stable  node, 

the  yawing  motion  can  approach  a  limit  cycle  of  circular  yawing  motion 

2 

or  diverge  to  large  values  of  K±  .  This  limit  motion  is  a  second 
property  of  non-linear  equations. 


* 


Although  Fig.  1  was  drawn  for  the  linear  equation,  the  amplitude  planes 
for  non-linear  equations  whose  only  singularity  is  located  at  the  origin 
are  quite  similar-  to  those  shown. 


If  neither  line  goes  through  the  first  quadrant,  the  amplitude  plane  is 
similar  to  that  shown  in  Fig.  lb. 
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If  the  zero  damping  lines  intersect  in  the  first  quadrant  and  0  ^ 

the  intersection  must  "be  either  a  node  or  a  spiral.  For  this  case  the  three 
singularities  on  the  axes  are  saddles  as  shown  in  Fig.  (3b).  When  the 
intersection  is  a  node,  it  lies  close  to  either  axis  and  so  this  case  is 
not  shown.  The  limit  cycle  predicted  for  a  stable  node  or  spiral  is  now 
a  limit  epicycle.  In  other  words,  for  a  wide  range  of  initial  conditions 
the  motion  should  eventually  be  epicyclic  with  certain  fixed  values  of 
amplitudes  and  no  damping.  Because  of  symmetry  considerations,  tpecspfecial 
case  of  a^Q  =01has  to  be  a  center  and  this  is  shown  in  Fig.  4. 

6.  COMPARISON  WITH  NUMERICAL  INTEGRATIONS 

After  developing  the  technique  for  treating  non-linear  fourth  order 
systems  of  a  certain  type,  it  is  natural  to  wonder  how  well  it  would 
predict  actual  motion.  Fortunately, the  exact  fourth  order  equations  of 
motion  for  a  missile  acted  on  by  non-linear  forces  and  moments  have 
already  been  programmed  for  the  NORC  computer  at  the  Naval  Proving  Ground . 

In  Reference  0  there  are  described  a  number  of  calculations  which  were 
made  to  investigate  the  effect  of  initial  conditions  on  the  character  of 
the  yawing  motion.  Since  these  calculations  involved  rather  large 
yarying  yaw  of  repose  and  a  non-constant  v,  the  comparison  of  theory 
with  these  calculations  will  be  deferred  to  the  appendix. 

Dr.  Cohen  and  Mr.  Hubbard  of  NPG  offered  to  make  a  number  of  special 

runs.jas  a  critical  test  of  the  equivalent  linearization  theory7  .  It  was 

felt  that  the  prediction  of  limit  epicycles  was  an  essentially  new 
prediction  and  that  this  prediction  should  be  checked.  Two  cases  of 
ten  NORC  runs  each  were  considered.  For  both  cases  the  Magnus  moment 

was  assumed  to  be  a  quintic  function  of  S  and  all  other  forces  and 

moments  were  assumed  to  be  linear. 

.  T  =  TQ  +  t2  e>2  +  T^  &V  (80) 
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are 


In  Reference  6  it  is  shown  that  the  effective  values  of  B^ 
r  In.  k  2  P  k 

[8JelaKl  +6kiK2  +5V  <8l) 

"  *2k  *  6KlV  +  5Kl4-  (82) 

• '  •  “i  *  “io  +  “i2(Ki2  +  2k22)  +  +  6k/  k/W'Jk/)  (85) 

a2  =  a20  +  a22(K22  +  2Ki2)  +  a2k^K2  +  ^±*2  +  5K/)  W 

.  _  t4 

where  . 

Thus  the  loci  of  zero  damping  are  now  hyperbolas. 

The  two  sets  of  aerodynamic  coefficients  and  6pin  were  selected  so 
that  the  left  tranches  of  the  zero  damping  hyperbolas  fell  in  about  the 
same  location  as  the  zero  damping  lines  in  Figures  3a  and  3b.  The  other 
branches  were  far  enough  to  the  right  that  they  could  be  neglected  for 
small  yaw  conditions.  The  exact  conditions  of  these  two  cases  are 
given  in  Table  1  and  the  location  and  type  of  their  amplitude  plane 
singularities  are  given  in  Table  2.  (Although  the  large  yaw  singu¬ 
larities  are  listed  in  Table  2,  these  points  are  not  important  in  this 
discussion. ) 

The  initial  conditions  for  Case  1  were  selected  to  lie  to  the  left 
of  the  dashed  curve  in  Figure  3a.  It  was,  therefore,  expected  that  the 
motion  would  quickly  become  a  pure  precessional  mode  with  amplitude .. 08l8 . 

p 

In  Table  3, initial  values  of  are  tabulated  with  their  values*  after 
18,000  calibers  of  travel.  From  Table  3,  we  see  that  the  final  values  are 
close  to  the  predicted  point.  In  Figure  5,  the  final  values  at  18,000 
calibers  are  plotted  on  a  greatly  enlarged  scale  to  show  how  the  motion 
is  slowly  approaching  the  pure  mode. 


The  amplitude  of  the  modes  may  be  easily  calculated  from  pairs  of 

successive  stationary  values  of  5.  (B  =  K,  +  .  B  =  I K  -  K  I  1 

^  max  ^1  *2'  min  1  *2 
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TABLE  1 

Parametric  Values  for  NORC  Runs 


*  =  .0055 
io3h  =  .too 


Jg=  0 


Both  Cases 

103  TQ  *  -  1.30 

103  t2  =  156 
103  t4  =  -  985 


Case  1  Case  2 

105M  =  -  4.16  -  2.08 

- — -  =  -  448  calibers  590  calibers 

TABLE  2 

Location  and  Type  of  Singular  Points 
Case  1 


K1 

*2 

Kl2 

K22 

Tyre- 

0 

0 

0 

0 

saddle 

.119 

0 

.0141 

0 

saddle 

0 

.0818 

0 

.006686 

stable  node 

.381 

0 

.1449 

0 

saddle 

0 

.390 

0 

.1523 

saddle 

.257 

.162 

.0660 

.0261 

stable  spiral 

Case  2 

0 

0 

0 

0 

saddle 

.115 

0 

.0132 

0 

saddle 

0 

.088 

0 

•  0075 

saddle 

.0268 

.0784 

.00072 

.00615 

stable  spiral 

.382 

< 

0 

.1458 

0 

saddle 

0 

.389 

0 

.1514 

saddle 

.250 

.169 

.0625 

.0286 

stable  spiral 
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TABLE  3 


N0RC  Results 

for  Case  1 

Initial 

Final 

Runs 

2  2 

10 

ioV 

io2k12 

2  2 
40  Kg 

1 

.0392 

.5083 

.0069 

.6512 

2 

.0400 

..5069 

.0046 

.6561 

3 

.4096 

.1823 

.0031 

.6577 

4 

.0655 

1.2343 

.0079 

.6464 

5 

.0610 

.5520 

.0058 

.6529 

6 

.0132 

.4761 

.0021 

•  6577 

7 

.0004 

.4900 

.0002 

.6642 

8 

.2601 

.1109 

.0015 

.6610 

9 

.0071 

1.2144 

.0061 

.6529 

10 

.0036 

.5256 

.0010 

.6626 

From  Table  2  we 

see  that  the  small 

yaw  limit  cycle 

for  Case  2  is 

zero -damped  epicycle  with  amplitudes  K±  =  .0268,  Kg  =  .0784.  Five  pairs 
of  initial  conditions  were  used.  Although  both  members  of  a  pair  bad  the 
same  initial  modal  amplitudes,  the  first  member  was  initially  at  the 
maximum  yaw  while  the  second  was  at  minimum  yaw.  The  equivalent  linear 
theory  makes  no  distinction  between  these  initial  conditions  and  the 
NORC  computations  seemed  to  verify  this  characteristic.  The  actual  initial 
amplitudes  are  given  in  Table  4. 

TABLE  4 


Initial  Squared  Amplitudes  for  Case  2 


Runs 

Kl 

K2 

2  2 
10K1 

l-f 

11,  16 

.0268 

.0784 

.072 

.615 

12,  17 

.0141 

.0784 

.020 

.615 

13,  18 

.0700 

.0500 

.490 

.250 

14,  19 

.0100 

.1380 

.010 

1.904 

15,  20 

.0194 

.0825 

* 

0 

OJ 

OD 

.681 
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Once  again  computer  runs  of  l8 , 000  caliber  duration  (30-1/2  maxima*  Of 
yaw)  were  made.  In  all  cases  the  motion  became  an  epicycle  with  zero 
damping  and  the  amplitudes  of  both  modes  agreed  with  the  theory  to  four 
decimal  places.  Since  runs  11,  12^.  15,  16,  17;  and.  20  had  initial  ampli¬ 
tudes  near  the  limit  amplitudes,  they  attained  their  limit  values  rather 
rapidly, i.e.,  in  less  than  5,000  calibers.  The  remaining  runs  required 
the  full  distance  of  18,000  calibers. 

In  Figure  6,  the  complete  amplitude  plane  for  Case  2  is  presented. 

The  small  yaw  position  together  with  the  initial  points  for  runs  (13,  l8) 
and  (14,  19)  are  shown  in  Figure  7*  The  time  history  of  the  amplitudes 
for  runs  13  or  18  were  calculated  from  Eqs.  (49  -  50 )  by  the  Exterior 
Ballistics  Analog  Computer.  The  values  of  the  two  modal  amplitudes  between 
successive  maxima  and  minima  were  computed  from  the  exact  N0RC  calculations 
for  Run  13  and  are  compared  with  the  analog  computer  results  in  Figure  8. 

The  agreement  is  quite  good. 

As  a  final  check  of  the  theory  it  was  decided  to  verify  the  predicted 
location  of  the  separatrix  for  Case  2.  (This  is  the  dashed  curve  in  Fig.  7 
which  separates  the  large  yaw  trajectories  from  the  small  yaw  trajectories.) 
For  a  value  of  K,-,  values  of  K,  were  chosen  which  bracketed  the  predicted 

o  2  1  2 

K  of  the  separatrix.  This  was  done  for  two  values  of  K2  •  Al.though.cthe 
theory  predicts  that  the  results  are  independent  of  relative  phase  of  the 

modes,  runs-  were  made  for  modes  both  initially  in  phase  and  initially  out 

? 

of  phase.  The  values  of  which  were  selected  are  given  in  Table  5. 

Exact  six  degree  of  freedom  calculations  were  then  made  on  the  NORC 
for  these  initial  conditions.  The  results  for  3000-caliber  long  tra¬ 
jectories  fell  into  two  groups.  For  one  group  the  yawing  motion  decreased 
while  for  the  other  it  grew.  For  one  run  it  was  not  possible,  on  the 
basis  of  the  first  3000  calibers,  to  decide  to  which  group  it  belonged. 
These  three  possibilities  are  denoted  by  S  (stable),  U( unstable),  and 
N  (neutral)  in  Table  5.  As  can  be  seen  from  that  table,  the  agreement 
with  the  theory  is  excellent. 
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TABLE  5 


Trajectories  to  Locate 

Separatrix 

C\J 

w00 

In 

Phase 

Out  of 
Phase! 

.006 

.00800 

3 

S 

.006 

.00825 

S 

S 

.00 6 

.00850 

S 

S 

.006 

.00875 

S 

S 

.006 

. 00900 

N 

U 

.006 

.00925 

U 

u 

2 

Predicted  : 

.00895 

.012 

. 00450 

s 

s 

.012 

.00475 

s 

s 

.012 

.00500 

s 

s 

.012 

. 00525 

s 

s 

.012 

.00550 

s 

u 

.012 

•00575 

u 

u 

.012 

.00600 

u 

u 

2 

Predicted  K,  : 

.00555 
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7.  SUMMARY 


1.  A  subset  of  the  set  of  fourth  order  linear  systems  with  constant 
coefficients  has  been  identified  and  its  important  property  of  an  epicyclic 
solution  described. 

2.  An  approximate  solution  has  been  obtained  for  a  class  of  non¬ 
linear  equations  which  can  be  linearized  to  members  of  this  subsbt. 

3.  This  solution  was  applied  to  the  equations  of  yawing  motion  of  a 
symmetric  missile  and  the  very  useful  concept  of  an  amplitude  plane  intro  - 
duced. 

4.  The  predictions  of  the  approximate  theory  have  been  compared  with 
exact  results  and  excellent  agreement  has  been  obtained. 
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APPENDIX  A:  DISCUSSION  OF  THE  INHOMOGENEOUS  EQUATION 

In  this  appendix  the  stability  characteristics  of  the  inhomogeneous 

equation  will  be  considered.  Since  for  constant  spin  a  small  aerodynamic 
.  12  , 

asymmetry  has  a  form  similar  to  the  gravity  term  in  Eq.  (4l),  the  treat¬ 
ment  will  be  general  enough  to  cover  both  cases.  More  precisely,  we  will 
consider  the  following  differential  equation.* 

+  (H+Jg  -  iv)X*  -  (M  +  iv  T)X  =  JX£  e1*  (Al) 

i*  i(vp  +  tJ 

where  JX_  e  T  is  either  G  or  J  X  e 
e  e  e 

co^d 


G  =  v  (— ^)  cos  0 
u 

iV  ^  ^  "  k2  ^ 
e  e 

J 

Kjj  is  the  asymmetric  force  coefficient 

Ky  is  the  asymmetric  moment  coefficient 
e 

X£  is  asymmetry  angle 

*0  initial  orientation  of  the  asymmetric  force 
0  Is  the  inclination  of  the  trajectory 

A  is  axial  moment  of  inertia  and  B  is  transverse  moment  of  inertia 

The  particular  solution  for  the  linearized  form  of  Eq.  (Al)  Kith  a  constant 

•  •  . ... 

has  the  form  e  where  0^'  =  ljr *  and  is  a  real  constant.  (For  the 
yaw  of  repose  case  \(/r  1  =  Oj  for  aerodynamic  asymmetry  \|f*  =  v.) 


*  £  * 
The  small  terms  -j - 


X*  and  iJL'X  have  been  omitted  for  simplicity. 
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An  important  difference  between  Eqs.  (A5  -  A6)  and  their  linearized 
versions  is  the  fact  that  more  than  one  solution  of  the  non-linear 
expressions  is  possible.  In  usual  practice  it  is  difficult  to  construct 
an  important  example  of  a  multivalued  yaw  of  repose .  For  missiles 
possessing  small  aerodynamic  asymmetries,  a  multivalued  trim  angle  is 
quite  possible.  As  is  shown  in  Figure  9,  which  is  based  on  Eq.  (A 6) 
for  a  common  non-linear  static  moment  curve,  the  well  known  linear 
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resonance  curve  is  quite  distorted  and  for  certain  spins  three  values  of 

are  possible.*  This  non-linear  response  curve  possesses Lthelcharacteristic 
non-linear  property  of  Jumps  between  roots  where  the  curve  possesses  a 
vertical  tangent. 


In  order  to  treat  the  more  general  case  of  mixed  oscillations,  we  will 
make  use  of  the  direct  substitution  method  of  Reference  6.  If  Eq.  (Al)  is 
written  for  cubic  non-linearities  in  M,  T,  and  J, 


+  (H+J  -iv)X' 


-  _Mq  +  «2 

i(t'P  +  tQ) 
e 

e  +  K,  e  J 

5 


+  iv  (TQ  +  T2 


62) 


(A7) 

(A8) 


2—222 
5  =  W  =  K±  +  K2  +  K?  +  K±K2 


±(01-02) 

i(02 -0-,) 

i(0  -0  )  i(0  -0  ) 

e 

+  e 

+K2K3 

e  J  +  e  J 

+  K3Ki  L' 


i(05-01)  i(01-05) 


+  e 


where  K.  =  e 
J 


-ajP 


0,  =  0iO  +  0,’P 

J  <JU  J  J, 


(A9) 


>  j  =  1,  2 


K,  =  constant  and 

5 

03  =  0jo  -•  v'p- 


It  is  quite  likely  that  the  middle  value  of  K,  is  unstable.  This  is  the 
case,  for  the  single  degree  of  freedom  case.  fs«&4,Refs.  n  and  1J.  ) 
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If  Eqs.  (A8  -  A9)  are  substituted  in  Eq.  (A7)  and  terms  of  the  same 
frequency  are  collected, 
i  0,  [“  o 

-L  I  /  .  jfll  ir  ,  /n  .  t  »  a  rA  f  f  xi  -t«.  m  \ 


K1  el  1  ai  +  ^i*)2  +  (H  +  Jg  -  iv  j ( -  a1  +  10^)  -  (Mq  +  iv  TQ) 

[*%» ' v. 


i(1ir0  "  ^30^ 


i0p  p  _  _ 

^  e  (-  a2  +  i0£»)  +  (H+Jg  -  iv)(-  a2  +  i02»)  -(Mq  +  iv  TQ) 

r*L-*  r^e1(*°-V 


:>  P 


e2  2 


J2  [fle32e 


i(  ^0  * 


10,  r  2 

+  e  0  J  -  t'  +  it*  (H  +  Jg  -  iv)  -  (M0  +  iv  TQ) 
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E  =1^2  +  iv  T, 


KxK2  e 


2  i(202-01)  2  i(20-0) 

e  +  Kx  K2  e  12 

2  2  i(2M0 

+  K2K5  e  5  *  +  K2^K5  e  25 

♦  KjKl2  el(20l’%)  +  ^  ,1(iW 

*  W  (e1<^-05)  *  e1'02^-0!)  +  /<W2r 


+  J2X£  e 


l(+0  +  +'P) 


K2K3 


i( 0,-0o ) 


+  KxK5  e 


i(0r02)  1(02-0,)“ 

+  KjKg  (e  1  2  +  e  2  1  ) 


If  certain  isolated  sets  of  frequencies  are  avoided*,  the  error  term  in 
Eq.  (A10)  vill  contain  frequencies  which  are  different  from  any  of  the 
"primary"  frequencies  0  ’ ,  Making  our  usual  assumption  that  the  solution 
to  the  non-linear  equation  has  the  same  form  as  the  solution  to  the 
linearized  equation,  we  neglect  the  error  term  E  in  Eq.  (A10)  and  replace 
the  local  amplitudes,  K^,  in  the  coefficients  of  cubic  terms  by  their 
average  values. 

'  *  ("  aj  +  ^j')2  +  (H  +  Jg  -  iv  H-atj  +  10j')2  -  (Mq  +  iv  Tq) 

-  («2  +  iv  T2)  [s2] 


Jej  '  J2 


e3j 


e  50  =  0  (All) 


(J  =  1,2). 


K?  e 


10 
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- 

J0  + 

J2 

HeJJ, 

e 

lt0 

t'(f’-v)  +  M0  +  M^b2] 

e5  +  1 

[v(To4T2 

V 

ej)-+’  <H+V 

(A12) 


The  cases  which  are  to  be  avoided  include  stability  factor  equal  to  one 
for  which  0-^ '  =  02'  and  resonance  for  which  either  0^'  =  01*  or  0^>  =  0^ 
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Eqs.  (All  -  A12)  may  be  extended  to  handle  more  general  non- linearities  by 

2  T  2k"| 

the  use  of  higher  order  polynomials  in  b  .  In  Table  6,  values  of  6  I  , 

r  2kl  ^  J  eJ- 

and  |J5  which  were  calculated  by  an  algorithm  similar  to  that  used 

in  Ref.  6,  are  tabulated  for  k  =  1,  2,  3* 


can  be  obtained  from 
subscripts  while 


-2k 

L  Jel 


and 


<-2k 

_6  Je3l 


e32 


by  replacing  1  by  2  or  3  in  the 


e35 


is  the  symmetric  part*  of 


&2klr 


An  important  feature  of  this  three  frequency  problem  can  be  seen 
from  Eq.  (A9).  According  to  this  equation  the  non-linearities  are 
functions  of  three  difference  frequencies.  These  difference  frequencies 
are  linearly  dependent,  and,  therefore,  this  treatment  is  essentially  an 
average  over  two  frequencies  and  suffers  from  the  hazy  physical  meaning  of 
such  a  process. 

As  an  application  of  Eqs.  (All  -  A12)  we  will  consider  the  amplitude 
plane  for  a  quintic  Magnus  moment  as  modified  by  a  constant  yaw  of  repose. 
If  Eqs.  (All)  are  modified  for  a  quintic  Magnus  moment,  it  follows  from 
their  imaginary  parts  that 


where  a 


JO 


aJ  -  “jo  +  aJ2 


(H  V  T0 


12  [b%3  +  aj4  l*4} 


ej 


(A13) 


20, 


-  V 


vT„ 


Q '  = 

j2 


2V 


-  V 


-  v  Ti, 


a 


J4 


2V 


This  symmetric  part  is  identified  by  braces  in  Table  6. 
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TABLE  6 


-P- 

-p- 


PJel  'l*!2  +  \2  *  Kj2  -  4  (K/  +  Kj2)  -  K/  4  2(K/  + 
^Jel  =^<K12  +  l^2  +  Kj2)  4  2(K  1%2  4  IC,2]^  4 


4  2(E^2  4  K/  4  K-^)  (Ku:!  +  Kj?)  4  te^Kj2 


2  .  „  2  .  „  2W„  2 


=  (Kxc  +  K2"  +  K3"){Kl"  +  3K2^  +  3K^)  +  2(K12K22  +  3K22K52  +  K^K^) 


LS°Jel  “  1(K12  +  *2  +  K32)  +  6(K/  +  K22  +  K32)(K12K22  +  +  K^2)  +  12K± 


2  2  2 
*2*3 


+  3(K12  +  K^  +  K?2)  {K^  +  K32)  +  l8(K12  +  +  K^)K^K^  +  3(K22  +  K?2  f  K± 


2,„2.,,2v2„2 


L"Je31  K3Xe 


2  2, 


e31  =  2K3Xe(Kl  +  ®2  +  K3  } 


L~  _le31  ”  ^Xe 


(Ki2  +  ^  +  K32)(K12  +  3k/  +  k/)  +  2K22(K12  +  K32)  +  K12K52J 


PJel  '  ^  +  ^  *  25R2 

Mel  -  M  +  6KlV  +  +  66R2(K12  +  “/>  +  3  C 

P]e2  ‘  ^  +  ^  +  2  8R2  “d 

pje2  -  (Kj,4  +  +  Xj1*)  +  6  6h2  (k22  +  at/)  +  3  6R\ 


In  most  cases  the  effect  of  T  on  8^  is  small. 


From  Eq.  (A5), 


v 


gd  cos  6 


(A14) 


Eqs.  (A13  -  All*-)  can  now  be  used  to  Interpret  the  results  of  Ref.  8. 

The  aerodynamic  coefficients  which  were  used  in  that  reference  are  listed  in 
Table  ff,.  Although  a  number  of  spins  were  considered,  we  will  limit  our 
attention  to  v  =  .0055.  The  quintic  expansion  of  the  Magnus  moment  was 
good  up  to  ll*°  and  so  only  the  small  yaw  singularities  are  listed.* 


TABLE  T 


Parameters  Used  in  Reference  8 

105H  =  .928  105Tq  =  -  2.83 

105Mq  =  -4.10  io5t2  =  313 

io\  =  -6  io5t^  =  -1976 


.0055 


=  1*52  calibers 


-  =643  calibers 


-  =  ll*96  calibers 

%' 


Since  the  frequencies  are  not  important  to  this  discussion,  the  indicated 
cubic  nature .of  the  static  moment  will  not  affect  our  results. 


45 


Small  Yaw  Singularities  for 
Initial  Yaw  of  Repose  of  .016 


2  2 

«2  V  V 

0  0  0 

0  .0144  0 

.082  o  .0067 


Type 

saddle 
saddle 
stable  node 


The  actual  variation  of  the  yaw  of  repose  as  a  function  of  p  is  given 
in  Figure  10.  Since  the  summit  value  of  the  yaw  of  repose  is  about  five 
times  its  initial  value  of  .016,  its  square  magnitude  grows  to  twenty- five 
times  its  Initial  value.  In  Figure  11  the  small  yaw  amplitude  plane  for  a 
constant  yaw  of  repose  of  .016  is  shown.  The  corresponding  amplitude  plane 
for  summit  conditions,  yaw  of  repose  of  .100,  is  given  in  Big.  12.  Thus  it 
can  be  seen  that  the  growth  in  6^  has  the  effect  of  moving  pure  mode  singu¬ 
larities  towards  the  origin  and  past  it. 

This  variation  in  6^  means  that  G  in  Eq.  (Al)  is  a  function*  of  p.  If 
G  and  therefore  changes  slowly  during  an  appropriate  period  of  the  motion, 
it  is  reasonable  to  expect  that  our  technique  will  still  apply.  According 
to  Table  7  the  shortest  period  associated  with  any  of  the  three  difference 
frequencies  i£  4^2  calibers.  From  Figure  10  it  can  be  seen  that  for  the 
first  thousand  calibers  the  yaw  of  repose  varies  reasonably  slowly  In 
comparison  with  this  period.  Since  the  results  of  Ref.  8  Indicate  that  the 
motion  over  this  distance  Is  sufficient  to  determine  the  character  of  the 
yawing  motion,  it  seems  reasonable  to  compare  the  predictions  of  the 
equivalent- linear  theory  with  those  of  that  reference. 

In  the  N0RC  calculations  of  that  report,  certain  physical,  parameters 
were  varied  and  a  large  number  of  trajectories  were  calculated  for  a 
variety  of  Initial  conditions.  It  was  found  that  the  computed  yawing 
motion  either  became  quite  large  or  damped  out.  The  critical  initial 
conditions  which  separated  these  two  types  of  yawing  motion  were  determined. 


*  Under  the  assumptions  of  Reference  8,  v  and  J  are  also  functiohs  of  p. 
Their  variation,  however,  is  overwhelmed  by  s  the  variation  In  6  .  In 
any  event  the  remarks  made  about  6  will  apply  to  these  terms  as  well. 
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The  initial  modal  amplitudes  associated  with  six  of  these  critical  sets  of 
initial  conditions  are  listed  in  Table  8  and  plotted  in  Figure  11.  Accord¬ 
ing  to  the  theory  these  points  should  fall  near  the  separatrix  associated 
with  the  nutational  saddle.  The  qualitative  description  of  this  saddle 
point's  motion  to  the  left  also  explains  why  the  point  for  Case  581  lies 
farther  to  the  left  of  the  separatrix  than  the  other  points.  In  all  events, 
the  agreement  is  quite  good. 

TABLE  8 


NORC 

Run  No. 

*1 

K2 

k22 

Standard 

.103 

.060 

.0106 

.0036 

Case 

446 

.106 

.054 

.0113 

.0029 

479 

.103 

.064 

.0107 

.0041 

504 

.101 

.062 

.0103 

.0039 

531 

.096 

.073 

.0092 

.0053 

581 

.079 

.105 

.0062 

.0111 
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APPENDIX  B:  CALCULATION  OF  THE  AMPLITUDE 
PLANE  FOR  NON- POLYNOMIAL  NON-LINEARITIES 

Although  the  results  of  this  report  have  been  stated  in  terms  of  general 
non-linearities,  the  applications  have  been  for  polynomial  functions.  These 
functions  possess  the  handicap  of  growing  rapidly  for  large  angles.  Most 
non#- linear  moments  either  level  off  or  tend  to  zero  with  increasing  angle. 
Since  constant  moments  imply  moment  coefficients  which  are  proportional 
to  &  ,  and  moments  tending  to  zero  can  be  described  by  higher  negative 

powers  of  &,  effective  values  of  negative  powers  of  &  are  of  some  interest. 


The  definitions  of  the  effective  values  of  &m  are 


2  2  2 

where  5  =  +  2^1^  cos  0.  If  m  is  a  negative  even  integer,  the 

integrals  may  be  computed  in  closed  form  by  means  of  Pierce  formulas 
Nos.  304-306.  This  has  been  done  for  m  =  -  2,  -  4,  -  6,  -  8,  and  the 
results  for  listed  in  Table  9*  The  results  for  W  he 

obtained  by  interchanging  the  subscripts  1  and  2. 

The  important  case  of  |p  J  e^,  which  arises  from  a  constant  moment, 
requires  the  use  of  elliptic  integrals.  After  some  algebraic  manipulation 
we  have  the  result: 


where  E±  is  the  complete  elliptic  integral  of  the  first  hind  usually  denoted 
by  K  and  E2'is  the  complete  elliptic  integral  of  the  second  hind  usually 
denoted  by  E. 
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TABLE  9 


2  2 
+  K1K2 


4„  2 


2„  4 


5  Kx  K2  -  3  Kx  K2 


2 

Since  the  inverse  powers  of  6  can  not  be  used  to  describe  aerodynamic 
moments  for  small  angles,  the  actual  moment  curve  would  have  to  be  approxi¬ 
mated  by  at  least  two  analytic  segments.  For  example  the  Magnus  moment 

could  be  represented  by  a  cubic  function  of  5, for  0=6=5,  and  a  constant 
^  ^  a 

for  -  6  =  6^.  The  Magnus  moment  coefficient  for  this  case  would  be 
determined  by  the  following  conditions: 


S  =  KT 
1o 


(B5) 


^  ' 


(B6) 


where 


V  -■ 


-1 


-  V  +  Va  B. 


o 


As  a  consequence  of  this  two  segment  approximation  of  the  Magnus  moment 
curve,  three  different  types  of  yawing  motion  have  to  be  considered: 

(1)  6  lying  in  the  interval  (0,  6&)  for  which  Eq.  (B5)  applies  and  the 
effective  values  of  ,6  can  be  used  to  describe  the  motion; 

(2)  5  lying  in  the  interval  (5&,  6^)  for  which  Eq.  (b 6)  applies  and 
the  effective  values  of  '6-1'  are  used;  and 

(5)  &  lying  in  both  intervals  for  which  both  equations  apply  and  the 

motion  calculated  by  means  of  Eqs.  (51  -  52). 

The  boundaries  of  the  corresponding  regions  in  the  amplitude  plane  are 


min 


(B7) 


&  v  =  K,  +  K~  =  & 
max  12b 

These  regions  are  shown  in  Figure  15. 

It  should  be  noted  that  the  amplitude  plane  for  arbitrary  non-linear 
moments  can  be  computed  by  numerical  means  from  Eqs.  (44  -  45).  This 
integration  would  be  simpler  than  that  of  the  complete  fourth  order  non¬ 
linear  equation  and  the  resulting  amplitude  plane  would  display  simple 
relationships  between  initial  conditions  and  the  type  of  yawing  motion. 
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APPENDIX  C:  AN  IMPROVED  FIRST  APPROXIMATION 


An  important  feature  of  the  K-B  equivalent  linearization  method  which 
was  developed  in  the  report  proper  was  the  use  of  an  epicycle  with  zero 
damping  as  the  first  approximation.  This  was  defined  by  the  equation 


/V 

X 


(Gl) 


where  %  =  0.JQ  +  ^  ’p 

-  a2  +|/a22  + 

2 

(V  2  'v  ’v  A-' 

a2  +  ^1  ^  0  an(l  a2  b^  are  the  values  of  a2  and  b^  when  X  and  X* 

are  zero. 

This  selection  of  a  first  approximation  suffers  from  the  handicap  that  it 
can  not  be  near  the  actual  motion  for  amplitudes  which  make  the  non-linear 
a2  and  b^  very  different  from  their  linear  values.  In  particular,  if  the 
small  amplitude  values  of  a^  and  b,  should  not  satisfy  the  inequality 
a2  +  4b  ^  y  0,  small  amplitude  epicyclic  motion  would  be  impossible  and 
the  first  approximation  would  not  exist.  Yet  it  is  quite  possible  that 
for  certain  non-linearities  large  amplitude  epicyclic  motion  could  occur. 

For  this  reason  we  will  consider,  in  this  section,  a  better  choice 
for  the  first  approximation.  From  an  examination  of  Eqs.  (38)  and  (40)  it 
can  be  seen  that  these  equations  could  be  considerably  simplified  if  a2  and 
b^  had  been  selected  to  be  the  average  values  of  a2  and  b^.  More  precisely, 
they  would  be  defined  by  the  relations 


2« 


0 


51 


Under  "these  definitions  Eqs.  (38)  and  (ij-O)  would  reduce  "to 


This  use  of  an  improved  first  approximation  which  is  based  on  def ig¬ 
nitions  (C2)  and  (C5)  would  clearly  overcome  the  problem  of  non-existent 
small  amplitude  epicyclic  motion*.  The  remainder  of  the  analysis  would  be 
identical  with  that  of  the  text.  The  calculation  of  the  amplitude  planes 
would,  however,  be  made  more  complicated  by  the  need  to  consider  the  vari- 
ation  of  with  amplitude. 


This  definition  also  possesses  the  advantage  of  yielding  the  exact 
frequency  of  constant  amplitude  pure  mode  motion.  This  can  be  seen  from 
the  fact  that  if  damping  is  neglected,  Eq.  (25)  is  linear  for  pure  mode 
motion  and  the  linear  frequency  formula  for  the  pure  mode  is  correct. 
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APPENDIX  D:  AMPLITUDE  PLANES  FOR  MISSILES  WITH 
NON-LINEAR  DAMPING  MOMENTS  AND  ZERO  SPIN 


A  rather  interesting  special  case  of  non-linear  damping  mddient  is  that 

f 

cL 


for  zero  spin.  Since  ,  Eqs.  (44  -  4-5)  assume  the  simple  form*: 


V 

K, 


2it 

r 


1 


J 


H(52) 


K2 

1  -  —  cos  0 

L  Kl 


a0  =  (D1) 


2rt 


V 

K2 


1 


<*> 

ri.h 

COS  0 

'  0 

L  2 

_ 

10  =  -  ^(K-l2,  Kg2 )  (D2) 


If  H  is  assumed  to  be  represented  by  a  power  series  in  6  ,  Eqs.  (D1  -  D2) 
become 


el 


(1 


10 


cos  0)  10 


The  small  gravity  term  and  small  ’  term  have  been  neglected  for 
simplicity. 


(D3) 


(H4) 
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The  definitions  of  J^S23*]^  should  have  a  simple  relationship  with  the  jjs2^^ 

listed  in  Table  I  of  Reference  6.  If  the  sign  in  fron^  of  the  coefficient 
of  the  brackets  which  appear  in  that  table  is  changed  from  plus  to  minus, 
this  new  effective  yaw  can  be  computed.  For  example 

'  *1  +  *22  -  M 


I?4]el  -  <K12  +  K/>2  +  -  2(Ki  +  *2®)  [k/| 


Ip  2  2  4> 

=  V  +2KiV- V‘ 


From  symmetry. 


He,  - 


(D5) 


(D6) 


(D7) 


'He,  '  ^  -  hk- 


(D8) 


Note.  If  H  had  been  assumed  to  be  a  polynomial  function  of  |  \»|2  =  k'k'  , 
interesting  simplification  would  follow.  For  0^  =  -  Eq.  (34) 


an 
becomes 


x'\'  = 


(X'X* ) 


K1  +K2  -  2KiK2  COS  (0^02 ). 

2* 


el=^  JQ  (Kl2+K 


(D9) 

(DIO) 


22  -  2^  cos  0)k  (1  -  ^  cos  0)  d0. 


Replacing  0  by  0  +  n  we,  therefore,  see  that 


(x-x- )k 


ej 


=  ,02k  &2k 


ej 


(DU) 
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Returning  to-"  the  assumption  that  H  is  a  function  of  8  we  now  consider 

the  properties  of  amplitude  planes  associated  with  Eqs.  (D3  -  D4).  Since 

these  relations  are  symmetric  in  the  two  modes,  we  expect  all  amplitude 

2  2 

planes  to  he  symmetric  with  respect  to  the  line  ^  •  The  first 

consequence  of  this  symmetry  lies  in  the  fact  that  the  origin  must  be  a 
node . 


In  order  to  consider  other  singularities  on  the  axes  we  expand  H  in  a 
Taylor  expansion  about  the  amplitude  K2 . 

H  =  HtK2)  +  H»(K2)  [s2  -  K2  (D12) 

(D13) 

(DU) 

2 

If  the  singularity  is  on  the  one  axis  at(K  ,  0),  vanishes  and  from  Eiq. 
(Dl3)it  can  be  seen  that 


ai  "  2 

‘  Hdc2)  +  E’dc2)  [Kl2  -  4 

a2  =  2  " 

H(K2)  +  H*  (K2 )  [K/  -  K ^jT 

HOC2)  =  0.  o  (D15) 

Shifting  the  origin  to  the  singular  point  by  the  translation  x  =  K-^  -  K2 , 

p 

y  =  K  ' ,  the  differential  equation  for  the  amplitude  plane  becomes 


dy  _  -y  (y  -  K2) 

"  =  X  (X  +  /) 


dx 


(D16) 


The  usual  tests  show  that  the  singularity  must  be  a  saddle.  Therefore, all 
pure  mode  singularities  with  the  exception  of  the  origin  must  be  saddles. 


The  case  for  a  quadratic  E  is  quite  simple.  The  damping  curves  are 
given  by  Eqs.  (D13  -  Dl4)  for  K2  =  0,; .E^  H(0]  "  and  =  E’(.0  ).  If 
Hq  and  Hg  have  the  same  sign,  the  origin  is  the  only  singularity.  (Fig. 
14a).  If  they  are  of  unlike  sign,  three  more  singularities  appear  in  the 


H, 


first  quadrant.  The  two  Sure  mode  singularities  at  (-  ■=-  ,  0)  and 

*2 


Hn 

(0,  -  —■)  have  to  be  saddles  while  it  can  be  shown  that  the  singularity 
*2 


at  (■ 


o 

— )  Is  a  node.  (Fig.  l4b). 

*2 
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For  negative  HQ,  this  second  node  is  a  stable  node  and  the  limit  motion 

rur 

is  planar  yawing  motion  with  amplitude,  Kx  +  K2  =  2  |/“  jT*  For  Planar 

yawing  motion,  the  yaw  equation  reduces  to  a  real  equation  in  the  magnitude 
of  yaw,  6. 


5’1  +  (Hq  +  Hg  62)B»  -  M  5  =  0,  M  4  0 


(D17) 


Eq.  (D17)  is  the  well-known  van  der  Pol  Equation  and  our  predicted  amplitude 
agrees  with  the  amplitude  obtained  by  other  methods.* 

If  a  quartic  H  is  considered,  singularities  can  appear  off  the  axes  in 

three  different  ways.  In  Table  10  the  off-axis-singularities  for  these 

three  cases  are  listed  and  the  amplitude  planes  for  cases  B  and  C  are  given 

in  Figs.  15  and  1 6.  (Since  the  amplitude  plane  for  case  Ai  is  essentially 

a  distortion  of  Fig.  lkb,  it  is  not  shown).  The  limit  motion  for  case  C 

is  particularly  interesting.  It  is  an  ellipse  with  semi -major  axis  r.  +  r, 

T'  3 

and  semi -minor  axis  r,  -  r.~ . 

4  5 


*  The  interesting  character  of  the  amplitude  plane  for  no  spin,  quadratic 
H,  and  4  0  was  pointed  out  to  the  author  by  H.  L.  Reed. 
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TABLE  10 


OFF-AXIS  SINGULARITIES  FOR  H  =  HQ 


\  5 


H. 


■2  '  2 


Case 

A 

Ho 

0 

-  4o 

\ 

2 

2  v 

2  >  r; 

.  ) 

node 

Case 

B 

1 

Pfel 

2  w 

.  h 

■  B 

W 

2> 

node 

;2> 

node 

Case 

_C 

<  0 


1 

*  Eo  ,  1 

V 

B 

_VJ 

<  h4  4  X 

-hm 

*2 


0 

Ek 


(r^2,  r^2)  center* 


center* 


*  According  to  the  criteria  of  Ref.  14,  this  point  is  a  center  if  the  cubic 
terms  which  appear  in  the  numerator  and  denominator  do  not  affect  the 
character  of  the  singularity.  Fig.  1 6  is  based  on  this  assumption. 
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FIG.  2 
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AMPLITUDE  PLANE  FOR  CUBIC  MAGNUS  MOMENT  (a|0 


6] 


10*  K  * 
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RUN  13  OF  CASE  2 
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